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Nonspherical stellar collapse to a black hole is one of the most promising gravitational wave 
sources for gravitational wave detectors. We numerically study gravitational waves from a slightly 
nonspherical stellar collapse to a black hole in linearized Einstein theory. We adopt a spherically 
collapsing star as the zeroth-order solution and gravitational waves are computed using perturba- 
tion theory on the spherical background. In this paper, we focus on the perturbation of odd-parity 
modes. Using the polytropic equations of state with polytropic indices n p = 1 and 3, we qualita- 
tively study gravitational waves emitted during the collapse of neutron stars and supermassive stars 
to black holes from a marginally stable equilibrium configuration. Since the matter perturbation 
profiles can be chosen arbitrarily, we provide a few types for them. For n p = 1, the gravitational 
waveforms are mainly characterized by a black hole quasinormal mode ringing, irrespective of per- 
turbation profiles given initially. However, for n p — 3, the waveforms depend strongly on the initial 
perturbation profiles. In other words, the gravitational waveforms strongly depend on the stellar 

)~~> , configuration and, in turn, on the ad hoc choice of the functional form of the perturbation in the 

\q • case of supermassive stars. 

| PACS numbers: 04.30.Db, 04.25. Dm, 04.30.Nk, 04.70. Bw 
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IT) ! I. INTRODUCTION 

O ■ 

Detection of gravitational waves is one of the greatest challenges in experimental and theoretical physics in this 
' decade. Several kilometer-size laser interferometers, such as TAMA [lj, the Laser Interferometric Gravitational Wave 
Observatory (LIGO) 0, and GEO are in operation now and VIRGO will be in several years. In addition to 
these ground-based detectors, the Laser Interferometer Space Antenna (LISA) with an arm length of 5 x 10 6 km has 
been proposed 00, and is planned to start taking observations in 2012. 

Nonspherical stellar collapse is one of the most promising sources of gravitational waves for both ground-based 
(3jT)[ detectors and space antennas. Ground-based interferometric detectors have a good sensitivity in the frequency range 
between ~ 10 and 1 kHz. Thus, the stellar core collapse of a massive star to a neutron star or a black hole is 
. — ' one °f the targets. According to 0, at the formation of a massive black hole, quasinormal modes of a black hole are 
excited and gravitational waves of high amplitude associated with the quasinormal modes are emitted. (See for 
a review of black hole quasinormal modes.) The frequency of gravitational waves associated with the fundamental 
5^ i quadrupole quasinormal modes of rotating black holes is 10] 

M x 



/ ~ (0.03 -0.08)^^(300- 800) (^—j Hz, (1.1) 

where the frequency is higher for more rapidly rotating black holes. Equation (|l.l|l indicates that formation of black 
holes of mass > 20M Q may be a promising source for laser-interferometric detectors. 

The frequency band of space antennas is between ~ 10~ 4 and ~ 0.1 Hz 0. This suggests that the formation 
of supermassive black holes may be one of the most promising sources. Although the actual scenarios by which 
supermassive black holes form are still uncertain, viable stellar dynamical and hydrodynamical routes leading to 
the formation of supermassive black holes have been proposed 0,0,0- In typical hydrodynamical scenarios, a 
supermassive gas cloud is built up from multiple collisions of stars or small gas clouds in stellar clusters to form 
a supermassive star. Supermassive stars ultimately collapse to black holes following quasistationary cooling and 
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contraction to the onset of radial instability fl4L ITa . Such dynamical formation of supermassive black holes may be 

a strong gravitational wave source for lisa m m E [HI23. 

The most hopeful approach for the computation of gravitational waves from stellar collapse to a black hole is to 
numerically solve the fully nonlinear coupled equations of Einstein and the general relativistic hydrodynamic equations. 
There has been much progress in this field in the last few years |2lj . However, it is not technically easy to compute 
gravitational waves with high precision in numerical relativity, since the amplitude of gravitational waves associated 
with the stellar collapse is not very large and as a result, they could be contaminated by numerical noises and/or 
gauge modes. To cross-check the numerical results and also to compute precise gravitational waveforms, it is desirable 
to have another method. 

As an alternative approach, linear perturbation theory has been developed [52, |U 01- I n this approach, wc 
decompose the fully nonlinear metric and matter field of slightly nonsphcrical profiles into a spherically symmetric 
dynamical field and linearized nonsphcrical perturbations. Because of progress in numerical techniques and compu- 
tational resources, Einstein's equations in spherically symmetric spacetime can be now accurately computed at low 
computational cost. Furthermore, due to the spherical symmetry of the background spacetime, the perturbations can 
be expanded in spherical harmonics. Thus all the equations for gravitational waves reduce to two simple (1+1) wave 
equations of even and odd parities, which can be numerically solved with high precision. Although this method is ap- 
plicable only to slightly nonspherical problems, the gravitational waveforms of small amplitude associated with stellar 
collapse can be computed with high accuracy. The result is also useful for calibration of fully nonlinear numerical 
results. 

The history of progress in linear perturbation theory of a spherically symmetric spacetime is as follows. The first 
work in this field was done by Regge and Wheeler j2f| and Zerilli [2(| . They derived the linear perturbation equations 
of odd parity |25j and of even parity pr| in Schwarzschild background spacetime. Subsequently, a gauge-invariant 
formalism of linear perturbations was developed by Moncrief |27| . Extending his work, Cunningham, Price, and Mon- 
crief |22j | derived perturbation equations on the Oppenheimer-Snyder solution for a collapsing uniform dust ball |28| . 
and computed gravitational waves emitted during the gravitational collapse of a dust ball to a black hole. Gerlach and 
Sengupta |29l | subsequently developed a gauge-invariant formulation of the linear perturbation on general spherically 
symmetric spacetimes. Using the Gerlach-Sengupta formalism, Seidel and co-workers 23] computed gravitational ra- 
diation from stellar core collapses, focusing mainly on the waveforms associated with the formation of neutron stars. 
They numerically solved the spherically symmetric general relativistic hydrodynamic equations using the May- White 
scheme [30L l3l[ . Harada et al. |3^ | studied scalar gravitational radiation from a collapsing homogeneous dust ball in 
scalar-tensor theories of gravity, using a similar method to that of Cunningham, Price, and Moncrief J^]. Iguchi, 
Nakao, and Harada [3^, 01 studied nonspherical perturbations of a collapsing inhomogeneous dust ball, which is 
described by the Lemaitre-Tolman-Bondi solution |3q]. Recently, Gundlach and Martm-Garcfa [24| have developed 
a covariant gauge-invariant formulation of nonsphcrical perturbations on spherically symmetric spacetimes with a 
perfect fluid, and derived coordinate-independent matching conditions for perturbations at the stellar surface. 

The dynamics of spherically symmetric spacetimes have been often studied using a method developed by May and 
White [30| . In this sche me, spherically symmetric spacetimes are described in terms of the so-called Misner and 
Sharp coordinate system |3q . in which a spacelike comoving slicing and an orthogonal time coordinate are adopted. 
This formulation is robust for the simulation of oscillating spherical stars and stellar core collapse to neutron stars. 
However, it is not robust enough to carry out simulations for black hole formation, because the computation often 
crashes before all of the matter is swallowed into a black hole due to inappropriate choice of the slicing condition [54j . 

To compute black hole formation, a null formulation proposed by Hernandez and Misner is well suited [3J]. In 
this formulation, spacetime is foliated by an outgoing null coordinate and thus the whole region outside the black 
hole horizon can be covered. The singularity avoidance of the null foliation is assured until the foliation reaches an 
event horizon if the cosmic censorship holds. Using this formulation, Miller and Motta |38l ] performed the numerical 
simulation of collapse to black hole formation. Baumgarte, Shapiro, and Teukolsky |3^| used this formulation to study 
neutrino emission in the delayed collapse of hot neutron stars to black holes. This formulation was also applied to the 
study of cosmic censorship [33. |40| and the formation of primordial black holes ^l| . Linke et al. 01 a l so computed 
the spherical collapse of supermassive stars using an outgoing null coordinate but different radial coordinate, i.e., 
Bondi metric to study the neutrino emissivity during the collapse. 

In this paper, we present a new implementation in linearized Einstein theory, in which spherical background 
spacetimes are computed with the Hernandez-Misner scheme, while nonspherical linear perturbations are treated 
using the single-null coordinate system. With the Hernandez-Misner scheme, it is possible to follow spherical stellar 
collapse to a black hole until almost all the matter has collapsed below the event horizon. The null coordinate system is 
well-suited for computation of gravitational waves emitted near the event horizon, which we want to study here. As a 
related work to the present treatment, Siebel et al. [43j presented simulations of gravitational collapse of neutron stars 
to black holes and the computation of quasinormal ringing in the spherically symmetric Einstcin-fluid-Klcin-Gordon 
system using an outgoing null coordinate without linear approximation. 
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This paper is organized as follows. In Sec. II we briefly review the evolution equations for spherically symmetric 
spacetimcs in terms of the Hcrnandez-Misncr null formulation. In Sec. Ill we describe the evolution equations for 
odd-parity gauge-invariant perturbations in the single-null coordinate system and then derive the explicit matching 
conditions at the stellar surface. In Sec. IV we explain the method for computation of gravitational waves in our 
gauge-invariant formalism. In Sec. V we describe numerical techniques adopted in the current implementation. In 
Sees. VI and VII, we present the numerical results of test simulations, and gravitational waveforms from the collapse 
of a supermassive star and neutron star to black holes, respectively. Section VIII is devoted to a summary. We adopt 
geometrical units in which G = c= 1, where G and c denote the gravitational constant and speed of light, respectively. 



II. BACKGROUND SPACETIME 



as 



A. 2 + 2 split of spherically symmetric spacetimes 

We decompose a spherically symmetric spacetime M into a product as M. = M 2 x S 2 . Namely a metric is written 

g^ u = diag(gAB, R 2 Jab), (2.1) 

where gAB, R, and ~/ a b are the (1+1) Lorentz metric, scalar function on Ai 2 , and unit curvature metric on S 2 . The 
greek indices /i, v, . . ., capital latin indices A,B,.. ., and small Latin indices a, b, . . . denote the spacetime components, 
M 2 components and S 2 components, respectively. The covariant derivatives on A4, M 2 , and S 2 are denoted as 
1^4 and :a , i.e., we define them from the conditions g^.x = 0, gAB\c = 0, and jab-.c = 0. 
The stress-energy tensor for general spherically symmetric spacetimes is given by 

V = diag(t A iJ, (t a a /2)R 2 lab ). (2.2) 

The totally antisymmetric covariant unit tensors e^B on M 2 and e ab on S 2 are defined as 

^ce CB = -ff/, ta C e cb = la b . (2.3) 

B. Hernandez-Misner formulation of general relativistic hydrodynamics 

We choose the Hernandez-Misner coordinate system of the form 

ds 2 = -e 2 ^du 2 - 2e^ +x ' 2 dudx + R 2 (d9 2 + sin 2 9d<p 2 ), (2.4) 

where x is a comoving coordinate, and ip 7 A, and R are functions of u and x. 

We assume that the stars are composed of a perfect fluid, for which the energy-momentum tensor is written as 

V = ( e + P) u n u v + P9^, (2-5) 

where e, p 7 and are the energy density, pressure and four velocity of the fluid. 
We define the following new variables 

U = e-^R, u , (2.6) 
T ee e-V*R, x -U= ^1-^ + U 2 , (2.7) 
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where m is the Misner-Sharp quasilocal mass 36]. Then, the field equations are written in the form |37 



U M = — 



1-ci 



P,x + 



e + p 



R? 



-A/2 



,-A/2 



(n) , u 



l-cl 

-e^4nR 2 pU, 

T + U 

R,x 
e -A/2 



U x 



2ur\ 



47Ti? 2J 
f 

n 



4:TTR 2 e x/2 (eT-pU), 

m + AirpR 3 



1 



TR 2 



(2.8) 

(2.9) 
(2-10) 

(2.H) 

(2-12) 

(2-13) 

(2-14) 
(2-15) 



n is the baryon rest-mass density, / = f(x) is an arbitrary function associated with the rescaling of the radial 
coordinate, and c s is the sound speed which is defined by 



„2 — 



dp 

7k 



s=const 



where s denotes the entropy. 

The regularity condition at x = gives the boundary conditions as 

R = 0, 
U = 0, 

r = i, 

m = 0. 

The boundary condition at the stellar surface x = x s is given by 

p = 0. 



(2-16) 



(2-17) 
(2.18) 
(2-19) 
(2.20) 



(2.21) 



In the original form of the Hernandez- Misner formalism [3jj, / is chosen to be unity. In this case, x coincides with 
the conserved mass /i contained in the interior to a shell. Another candidate for / is 



fix) 



4irn(uo, x)x 2 

(T+u)(u Q , x y 



(2.22) 



In this case, x coincides with the circumferential radius of the shell on the null surface u — uq. In this paper we adopt 
Eq. H2.22|l for /, since it has a nicer feature for the integration of nonspherical perturbations which we will describe 
in the next section. 

We assume that the exterior of the star is a vacuum. Then, the zeroth-order solution is the Schwarzschild spacetimc 

as 



ds z 



2M\ 



dT 2 



2M\ 1 



dR 2 + R 2 (d6 2 +sin 2 



(2.23) 



where M is the gravitational mass of the system. To compute gravitational waves in this background, it is convenient 
to introduce null coordinates u and v defined as 



v. 
v 



T-R*, 
T + R„ 



(2.24) 
(2.25) 
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where 



Then, the line clement is of the form 



R, = R + 2Mlnl — -1 ] . (2.26) 



ds 2 = - (l - dudv + R 2 (d9 2 + sin 2 6W0 2 ). (2.27) 

In the following, we refer to the standard outgoing null coordinate u as an observer time coordinate according to |39| . 

The ratio of the proper time interval dr s of an observer at the stellar surface to the observer time interval du is 
given by 

= (T + U)(u,x s ). (2.28) 

du 

The lapse function a at a; is defined by the ratio of the proper time interval dr of each fluid element to the observer 
time interval du as 

. dr dr dr s e^ u ^ ,„ . . 

a {u,x ) = — = — — ? = — 7 T (T + U)(u,x s ). 2.29 

y ' du dr s du e^("> x <>) y ' 

We note that the lapse function a is directly related to the observed redshift z as a = 1/(1 + z). 

In solving the dynamics of a spherical star, the boundary condition for -0 is arbitrary. In the present computation, 
we choose e^(u, x s ) = 1 for the boundary condition of ip. Then, we can identify the null coordinate u with the proper 
time t s of the comoving observer at the stellar surface. 

III. ODD-PARITY PERTURBATIONS 
A. Gauge-invariant perturbations 

Perturbed metric and matter fields of odd parity are denoted by 



A ' 9 ^ " I hAS a h(S^+ a S b:a ) ) ' (3 ' 1} 

= ( Ai.Q a+(q , i c, ^ )' (3-2) 



At A S a At{S a:b + S b:a ) 

where Y, S a = e a b Y- b and S a -.b + Sb- a are the scalar, vector and tensor harmonics, respectively. Here, the suffices I 
and m are omitted for simplicity. The scalar harmonic function Y satisfies 

l ab Y :ab = -1(1 + 1)Y. (3.3) 

The gauge-invariant perturbations are defined as 

kA = h,A — h\A + 2hvA, (3.4) 

L A = At A -Qh A , (3.5) 

L = At-Qh, (3.6) 

where 

va - (3-7) 

L is identically zero for 1 = 1, and no perturbation of odd parity appears for I = 0. 

For a perfect fluid, the matter perturbations of odd parity are only specified by the four- velocity perturbations as 

A Ufl = (0,f3S a ), (3.8) 

where (3 is a function of a; and u and completely determines the matter perturbations. The concrete form is determined 
by solving the field equations [cf. Eq. I|3.10|l ]. In terms of /3, the gauge-invariant perturbations are written as 

L A = /3(e + P )u A , (3.9) 

and L = 0. 
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B. Field equations: interior 



The covariant field equations for nonspherical perturbations in the stellar interior were derived by Gerlach and 
Sengupta [2j| for general matter fields and have been recently reformulated by Gundlach and Martrn-Garcia [24| for 
a perfect fluid. We follow |24| to derive the basic equations for the perturbations. 

The perturbation equation for the matter field, A (#*„.„) = 0, is integrated to give 



= 



-A/2 



R 2 (e+pY 



where j = j(x) is an arbitrary function of x. Integration of j by dx yields a conserved quantity J as 



J 



jdx, 



where § JS^ sm0d0d(j) corresponds to the z component of angular momentum for 1 = 1. 
At R = 0, (3 should satisfy the regularity condition as 

(3 = R l+1 fi, 

where j3 is a regular function at R — 0. The regularity of (3 also leads to the following condition 

e -x /2j = R i +3]j 



(3.10) 



(3.11) 



(3.12) 



(3.13) 



where j is a regular function at R = 0. In terms of the function j, the gauge-invariant matter perturbations are 
described as 



Thus, in the following, we specify j instead of (3. 

The metric perturbations of odd parity are characterized by a master variable as 



= e -i>-\/i 



n = e 



k r 
R? 



R? 



(3.14) 
(3.15) 

(3.16) 



where II is related to a variable tti introduced in j2^| as tti = 1(1 + 1)R 2 H. II should satisfy the regularity condition 
at R — as 



n = r'- 2 u, 

where II is a regular function at R = 0. n satisfies the following wave equation for / > 2 



(3.17) 



2(1 + 1) 
R 

-(1 + 2) 



Me - p) + (I - 2) 



2m 



n 



16tt 



'Re-4>-W( e *j) + {(l + i)r + 2U}j 



The relation between Ua and II is given for I > 1 by 

(l-l)(l + 2)k u = 16^i? 2 L u + (i? 4 n), u -e^- A / 2 (i? 4 n)^, 
(l-l)(l + 2)k x = IQitR 2 L x - (R 4 U), X . 

For I > 2, the gauge- invariant metric perturbations are obtained from Eqs. (|3.18|l . 13.19|l and H3.2t)|l . However, for 
1 = 1, Eqs. jSJUl and (|3~^U|) give 



(3.18) 

(3.19) 
(3.20) 



i? 4 n = 16?r J, 



(3.21) 



where we assume that the perturbation is regular at R = 0. Equation l|3.21|l implies that there is no gravitational-wave 
mode for 1 = 1. 



7 



C. Field equations: exterior 



We have two metric perturbations kx and k R in the exterior. Defining the master variable II as 

n=(¥) , (3.22) 



\R 2 J ,t V-^ 2 / ,r 

we find the following wave equations for I > 2 

_$ TT + $ ^ _ V (R)$ = 0, (3.23) 

where 

$ = i? 3 n, (3.24) 

nR)s ( 1 _-)(Mi±ll_-). 

Here, we note that $ is related to the variable ip defined in as ij) — 1(1 + 1)$. Using the double- null coordinates 
(u, v), we obtain 

4$, os + V(R)& = 0. (3.26) 
Equation (|3.23|) has a static solution with an appropriate fall-off at infinity as 

q (2M\ l ( 2M\ 

where F(a, b, c; z) denotes the hypergeometric function and q, which has the dimension of length, corresponds to the 
multipole moment of the system. Because [1(1 + is factorized out in the above equation, the definition of the 

moment is the same as that defined in 22]. This static solution is used for providing initial conditions of metric 
perturbations (see Sec. V). 

The relation between kA and II is given for I > 1 by 

(i-\)(i + 2)k T = -(i-^f) (r*ti) iR 



R , 

(i? 4 n), ffl - (i? 4 n) iC , (3.28) 



(l-l)(l + 2)k R = -(!--_) (J2*n) 

2M X 



2M N 1 



,T 



= -(1-^) [(RiTl), a + (R*Tl)M- (3.29) 



For I = 1, we find the solutions of Eqs. I|3.28|l and 

i? 4 II = const. (3.30) 

The integration constant in Eq. I|3.30|l is related to the total angular momentum of the fluid perturbation for 1 = 1. 
The master variable II is time independent in the exterior for I = 1 as shown in Eq. H3.30|) . This constancy implies 
the conservation of the angular momentum of linear perturbation in the spherically symmetric background. 

To compute nonspherical metric perturbations, we divide the background spacetime into regions I— III (see Fig. [T]). 
Region I is defined as the interior of the star. Region II is an intermediate exterior region from which one can emit 
an ingoing null ray that encounters the stellar surface before the stellar surface is swallowed into the event horizon. 
Region III is the exterior region outside region II. Region II is introduced to help the matching procedure at the 
stellar surface in numerical computation. Hereafter, we will refer to the ingoing null surface which divides the exterior 
regions into regions II and III as the junction null surface. Such an elaborate procedure is needed to calculate the 
late-time gravitational radiation extracted at the point far from the star and to assure that the matching condition is 
satisfied at the stellar surface simultaneously. 

For the convenience of computation, we introduce new null coordinates u and v in region II. We identify these null 
coordinates u and v with the values of the proper time r s of an observer comoving at the stellar surface x — Xg cit its 
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intersection of an outgoing ray (u = const) and an ingoing ray (v = const), respectively. Namely, the stellar surface 
is given by u — v. 

When we define functions A and B as 

A(u) = d i = ¥ ^j(u, Xa ), (3.31) 
B(v) = ± = f ± l! (v,x t ), (3.32) 

we can rewrite the wave equation for $ as 

4$ sc + A(u)B(v)V(R)<& = 0. (3.33) 

We integrate Eq. I|3.33|l in region II. The event horizon is given by a finite value of u. It is found that the effective 
potential term V = A(u)B(v)V(R) is regular on the event horizon in this coordinate system, which helps numerical 
integration of the wave equation. 

D. Matching 



The matching condition at the stellar surface x — x s for the odd-parity perturbations is obtained from the continuity 

A 

6abu B .24]. The explicit equations are 



condition for II, u a La, and n^IIu — 16itR 2 u a La for I > 1, and for n A kA and u a Ua for I > 2, where ua 



n in = n cx , (3.34) 

-e-^II in , u + e- A / 2 n in , x - 167ri?- 4 e- A / 2 j 

"Ilex g + II ex y. (3.35) 



T+U T-U 

Using u^IIjui^ = u A H on t\A, we can derive an alternative form of the matching condition as 

-2e-^n in ,„ + e- A / 2 n ia , x - 16nR' 4 e-^ 2 j 



= = -2n cx , s . (3.36) 
For 1 = 1, the matching conditions lead to 

n in = IStt^, (3.37) 

n cx = 16tt^. (3.38) 

Thus, the gauge-invariant variable II is completely determined by the initial distribution of perturbed angular mo- 
mentum in the star. 

IV. GRAVITATIONAL WAVES 

To compute gravitational waves in the wave zone, it is convenient to adopt the radiation gauge. In this gauge, the 
following tetrad components denote the + and x modes of gravitational waves; 

h + = \{h 66 -h u ), (4.1) 

hy = h H . (4.2) 
Hereafter, we adopt the following choice for the orthonormal bases: 

Yio for m ~ and 

Y hm± ee -j=(Yi, m ± Y h _ m ) for m + 0. (4.3) 
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The metric perturbations h + and h x in the radiation gauge are written in terms of as 

h+ = + 2) ^ + co™t)A+(9, <f) + 0{R- 2 ), (4.4) 

hx = ^ ~ 2) ^(<& + const)^ (0, 0) + 0(R- 2 ), (4.5) 

where the angular dependence can be explicitly calculated as 

■A+(0,(f>) = Se-.e . 2 n S<t>-.<t> 

sm 9 

= 2 {^~a Y M ~ —7^^") . ( 4 - 6 ) 

AA6,ct>) 



sm# ' sin & tan ( 



sm( 



2 -Y e -2m 2 — i^-Y + ^ + ^Y (4.7) 



tan ' sin 9 



Here, Y denotes one of the bases shown in Eq. 14.311 . The luminosity of gravitational waves Pi. n for each mode in 
terms of the master variable is given by (see [22, 133] for m — and also |4j 



- ^m^? ^ Al+A)VA% (4 - 8) 



where the subscript n denotes and m±. 



V. NUMERICAL METHOD 
A. Numerical integration 



The spherically symmetric stellar collapse is computed using the single-null comoving coordinates. Our method 
is essentially the same as that used by Baumgarte, Shapiro, and Teukolsky jS^. An artificial viscosity term is 
incorporated to deal with shock waves. The details of numerical method which we adopt are found in . 

To solve the perturbation equations for the interior of the star (region I), it is convenient to decompose Eq. (|3.18l) 
into the first-order differential equations as 

~2e-^- x ' 2 Q ,„ + e -^- A /2( e ^-V2 Q) ^ + Hii±ll{_(r + U)e~*P + Te~ x/2 Q} 

R 



-(1 + 2) 



47r(e - p) + (I - 2) 



2m 



n 



16tt 



-2e 



Re-^- x / 2 (e^j), x + {(I + l)r + 2U}j 
P,x+e 



^ - /> , !( - A / 2 (e^~ A / 2 Q) x J- 2 ^ + jj 



-(1 + 2) 



4Tr(e-p) + (l-2) 



2m 



R 

n 



{-(r + Z7)e-^P + re- A / 2 Q} 



= 16tt 

n,„ = p, 
n >x = Q. 



Re-^- x/2 (e^j), 



R? 

{(l + l)T + 2U}j 



(5.1) 



(5.2) 

(5.3) 
(5.4) 



We note that the regularity at the center requires 



1>-\/2, 



(5.5) 
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Equation i|5.1[l constitutes a hyperbolic-type partially differential equation to which we apply the first-order up- wind 
scheme to stably evolve the function Q. Other equations constitute ordinarily differential equations, and thus the 
integration is carried out with the second-order Runge-Kutta method. The Courant-Friedrich-Lewy (CFL) condition 
for the stability of integration limits the n-th time step Au n as 

Au n = 2C*mine~^" +( ^ /2) Aa; 4 , (5.6) 

i 

where C(< 1) is the Courant number. 

In the exterior of the star (regions II and III), the double- null coordinates are adopted. In region II, Eq. H3.33[l is 
decomposed into the first-order differential equations as 

Z,v = ~A(u)B(v)V(R)$, (5.7) 

Wu = -jA(u)B(v)V(R)$, (5.8) 

S jfi = Z, (5.9) 

- W, (5.10) 

and in region III, Eq. (|3.2(i|l is decomposed as 

Z« = -\v(R)<S>, (5.11) 

W, u = ~V(R)$, (5.12) 

= Z, (5.13) 
<&,v = W. (5.14) 

For integration of these equations, we use the finite differencing scheme proposed by Hamade and Stewart [45j. 

The coordinates [u, v) depend on the spacetime trajectory of the stellar surface. Thus, they are determined after 
the spherically symmetric stellar dynamics is solved. For this reason, we divide a numerical simulation into two steps. 
In the first step, we carry out numerical computation for the zeroth-order background solution taking into account 
the CFL condition for the first-order nonspherical perturbations, and in the second step, we evolve the nonspherical 
perturbations. 

For accurate numerical integration, the distribution of grid points in region I plays a quite important role. We 
adopt the equally spaced grid in terms of the initial circumferential radius. In region II, the distribution of grid points 
is automatically determined in computing zeroth-order solution. In region III, the equally spaced grid in terms of u 
and v works well. 



B. Matching 

At the stellar surface, the matching conditions (|3.34l) and (|3.3ti|) for II and its derivatives determine the boundary 
conditions for S, P, and Q at x — x s in region I, and for Z and W at u = v in region II. In matching solutions of 
regions I and II, we have to be careful since the outermost several mass shells form a very disperse envelope during 
the collapse. This implies that the grid resolution around the outermost envelope is not very good. If the matching 
between regions I and II is done at the location of the outermost mass shell, a large numerical error is produced. 
To suppress the numerical error, we match the solutions at the location of a mass shell which is not outermost one 
but is located slightly inside the outermost one. We find that this works quite well. We have also checked that the 
numerical results do not depend on the location of the matching if more than ~ 95% of the total mass is contained 
inside the mass shell for the matching. 

At the junction null surface, the matching is not necessary, since the double-null coordinates are adopted in both 
regions. We only need to store the numerical data set on the junction null surface in region II, and use it as a part 
of initial data for the integration of region III. 



C. Initial data 



In the present formulation, the null cone composed of u — uq, u — uq, and u = uq should be taken as an initial null 
surface. Since it is the characteristic surface of the wave equation, we only need to specify IIi n it(i?) and /3j n it(-R) there. 
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We provide the initial data set in the following manner: In region I, we first give the functions IIi n it(-R) and Pmit(R)- 
Then, Tl(uo,x) and j{x) are specified from Eqs. (|3.10|) and 13.17|) . From the initial data TL(uq,x), we can determine 
Q(uo,x) on the initial slice by differentiation through Eq. (|5.4I) . Then, from fl(u ,x), Q(u ,x), and j(x), we can 
obtain P(uq, x) by integrating Eq. I|5.2I) with the central value given by the boundary condition (|5.5I) at the center. 

In region II we first provide $>(uo,v) = Tl(uo,v)R 3 (uo,v). Then, W(uq,v) is obtained by differentiation through 
Eq. (|5.1Q(I . and Z by integrating Eq. (|5.7|1 with the initial value given by the matching condition 1)3. 36|) at the surface. 

In region III the method for construction of the initial data sets of W, and Z is the same as that in region 
II, except for that the initial values needed for integration of Eq. (|5.11|) are given at the point on the junction null 
surface. 

When the background spacetime is initially momentarily static, it is natural to choose momentarily static initial 
data sets. Thus, the initial data set is given in the following manner. First we specify Pi a it(R), and compute j(x) 
from Eq. i|3.10[l . Then, we determine TL(uo,x) and Q(uo,x) by integrating Eqs. 1|5.2|) and 1)5. 4f> using the condition 
P{uq,x) — P, x (uq,x) — with the initial value II(uo,0) and Q(uq,0) = 0. Here, we have to tune II(uo,0) so that it 
matches the static exterior solution at the stellar surface. The matching is written as 

fl - 16n]R = 0, (5.15) 

where F(l - 1, 1 + 3, 21 + 2; 2M /R) is abbreviated as F t . 

After this procedure, the multipole moment q is determined by matching the interior solution with the exterior 
static solution l|3.27[l . Then, using the moment q, the initial data set for the exterior solution is provided by Eq. i|3.27[l . 



-A/2 



Q 



(21 + 1) + 



2M F( 
~R~F~i 



D. Event horizon 



The effective potential term of the wave equation for II is regular on the event horizon. However, it contains a 
term of the form "zero divided by zero" (for R — > 2M, A — > oo, and V — > 0). This implies that it is not possible 
to numerically integrate the wave equation for n on the null surfaces which are very close to the event horizon. To 
integrate the wave equations until the null surface reaches the event horizon, we use an extrapolation for the value of 
n on the junction null surface in the neighborhood of the event horizon as 

ttEH _ n 7V max 

11 = nAfmaX + R^_ RNm J R - RN ™)> ( 5 - 16 ) 

where the value n EH on the event horizon is estimated as 

n EH = ^ + |__E__ (jRE h _ RN ^ (5 17) 

pjiv max anc j jj>JV max are values for n and R on the junction null surface at the iV max th (final) time step in region II 
and i? EH = 2M. To compute gravitational waves in region III emitted in the neighborhood of the event horizon, this 
extrapolation is necessary and quite helpful. 



VI. CODE TESTS 



A. Spherically symmetric hydrodynamic code 

As the first test of the spherically symmetric general relativistic hydrodynamic code, we performed a simulation for 
the homogeneous dust collapse. The numerical solution can be compared with the Oppenheimer-Snyder solution [28| . 
We chose the total grid number as 500 in this test. As the initial condition, we set a momentarily static dust ball of 
uniform density and of radius AM . Precisely speaking, the dust surface is maximally expanded on the initial null slice. 
In this case, the event horizon is located at r s = 7.243Af, where r s denotes the proper time of an observer comoving 
with the stellar surface. 

In Fig. |21 we show the time evolution of the coordinate velocity profiles as a function of circumferential radius at 
selected time steps. The crosses denote the numerical results and the solid curves the exact solutions. The figure 
shows that the numerical results agree with the exact solutions within 0.01% error. We also note that the computation 
can be continued until the null hypersurface reaches a surface very close to the event horizon. 
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To demonstrate that the code works well even in the presence of pressure, we carried out a long-term evolution of 
spherically symmetric stars in equilibrium. Here, we adopt the polytropic equations of state as 

e = n(l + e), (6.1) 
p = Kn r % (6.2) 

ne = - P _ - , (6.3) 

where e is the specific internal energy, and T a an adiabatic constant, for which we choose as 2 and w 4/3. With 
T a = 2, moderately stiff equations of state for neutron stars are qualitatively approximated. With T a 4/3, the 
equation of state for supermassive stars of mass > 10 6 M Q , in which the radiation pressure dominates over the gas 
pressure, is well approximated. According to |15| . the equation of state for supermassive stars may be approximated 
by the polytropic equation of state with T a as 



r a = t +0.00142 (*L- ] «U) 



-1/2 

3 ' V 1O6M 0, 

Thus, we adopt T a = (4/3) + 0.00142 to model a supermassive star of mass 10 6 M Q . The equilibrium configurations 
are obtained by solving the Tolman-Oppenheimer-Volkoff equation. 

In Fig. |3 the gravitational mass as a function of the central density for the equilibrium configurations with T a = 
(4/3) + 0.00142 and 2 are shown. For simulations in this paper, we adopted four static configurations of nearly 
maximum mass along the equilibrium sequences (see A-D in Fig. |3J) as initial data sets. Models A and C are stable 
against gravitational collapse, while B and D are marginally stable. We summarize these models in Table [I] 

TABLE I: Models for equilibrium solutions. The unit ofc=G = M = lis adopted. 



Model 


r a 


Central energy density 


Radius 


A 


(4/3) + 0.00142 


1.852 x 10" 9 


1.900 x 10 3 


B 


(4/3) + 0.00142 


3.707 x 10~ 9 


1.508 x 10 3 


C 


2 


5.944 x 10~ 3 


5.501 


D 


2 


1.044 x 10" 2 


4.745 



For long-term test simulations, we adopt models A and C. To induce a small oscillation, we initially increased the 
specific internal energy e uniformly by 1% both for model A and for model C. In Fig. 0] we show the time evolution 
of the central density. It is found that the stars oscillate in a periodic manner. By performing the Fourier analysis, 
we measured the period of this stellar oscillation, and found that the oscillation angular frequency is 



M 1 ' 2 



u>- 



(6.5) 



i? 3 / 2 

where Q = 9.59 x 10 -2 for model A and 0.788 for model C, where the initial stellar radii are given by 1.900 x 10 3 Af 
and 5.501M for models A and C, respectively (see Table I). 

According to a post-Newtonian theory, the angular frequency of the radial oscillation for supermassive stars can be 
estimated as 

^ = — — (r a - r crit ) , (6.6) 

where W , I, and r cr i t are the gravitational energy, moment of inertia, and critical polytropic index, respectively. 
Here, r cr ;t is estimated in the post-Newtonian approximation as |46J 

4 . „_/2M\ 



r crit = - + i.i25^— y (6.7) 

From this analysis, ui should be 9.62 x 10 -2 for model A. Thus, the numerical result agrees with the analytic one 
within 0.5 % error. 

The angular frequency of the radial oscillation for neutron stars can be also calculated using the semianalytic formula 
derived by Chandrasekhar |46j . From this approximate formula, we can predict the oscillation frequency should be 
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~ 0.9/c»y 2 01 an( i hence u> ~ 0.063Af — 1 . On the other hand, the numerical result indicates that to ~ 0.0616A/ -1 , 
which is again in good agreement with the analytic result within a few percent error. 

Models B and D are marginally stable against gravitational collapse. Thus, if the internal energy is reduced, they 
start collapsing. We have checked that they indeed collapse to black holes. The detailed results of the gravitational 
collapse are described in Sec. IVIII 



We have carried out a wide variety of test simulations for our perturbation code. First, we checked that gravitational 
waves accurately propagate on the flat Minkowski background. For 1 = 2, the exact solutions for linear gravitational 
waves in the Minkowski spacetime are obtained as [48| 



Since the background spacetime is flat, the matching surface between regions I and II is artificial. As the matching 
surface, we chose the surface of r = 1 in this test. In Fig. [5] we show numerical results with the exact solution for the 
waveforms observed at r — 5. To demonstrate the convergence of the numerical solutions, we performed simulations 
with three levels of the grid resolution as Ar = Au = Av = 0.001, 0.002, and 0.004. Figure |3] shows that numerical 
results converge at first order to the exact solution. We note that the number of grid points per wavelength is ~ 1000 
for the best-resolved case. 

Next, we computed the propagation of gravitational waves for I = 2 in the collapse of a homogeneous dust ball. 
The same type of computation was already carried out , so that we can calibrate our code by comparing the 

present results with previous ones. 



In |22(, the momentarily static perturbations are provided at the initial hypersurface. Thus, we also choose the 
momentarily static initial data on the null hypersurface (i.e., the dust surface is assumed to reach the maximum 
expansion at initial slice). However, the word of caution is appropriate here. The previous computation was carried 
out using 1+1 coordinate system (not single- null coordinate system). On the other hand, we choose the single-null 
coordinate system. Namely, the coordinate system and also the time slicing are different between the two. Thus, 
even if we set the same function of /3 in different coordinate systems, this does not imply that we give the identical 
initial condition. Moreover, the former formulation has ambiguity in determining the second-order time derivative of 
gravitational perturbation, while the latter one does not. Hence, the momentarily static initial data sets for these two 
formulations are different from each other. It implies that it is difficult to precisely compare the results obtained in 
these two formulations. 

To calibrate the effect of different initial conditions on the results, we provided two perturbation profiles fixing q. In 
the first one, we gave /?i n it(-R) = const on the initial hypersurface, and in the second one, /3i n it(-R) °c exp[— (i?/i? c ) 2 ], 
where R c is chosen to be one-third of the initial surface radius. 

In Fig. the total radiated energy of gravitational waves are summarized. In these numerical computations, 
gravitational waves are extracted at R = 40Af. For comparison, we also plot the results of [22j. We note that in this 
figure, q is normalized so that q = 2M. It is found that numerical results of two different profiles of j3 differ by a 
factor of ~ 1.4, and that the results of are greater than our results by a factor of ~ 3. However, the total radiated 
energy systematically decreases as the initial dust radius increases for all the cases in the same manner. Thus, we 
conclude that our numerical results agree with the previous ones besides a possible systematic error associated with 
the difference of the coordinate conditions. 

The waveform, luminosity, and integrated total energy of gravitational waves from a collapsing dust ball with the 
initial radius R = 20 M for I = 2 are plotted in Fig. [7| Gravitational waves are extracted at R = 40M. As found 
in p^ . the quasinormal ringing oscillation and subsequent power-law tail characterize the gravitational waveforms. 
The complex frequency of the fundamental quasinormal mode was calculated to be 2Mw = 0.74734 + 0.17792i by 
solving the eigenvalue equations |49j| . On the other hand, our numerical results show that the complex frequency of 
the damped oscillation is given by 2Mu> = 0.752 + 0.179i. Thus, the numerical results agree with the theoretical value 
within 1 % error both for the frequency and for the damping rate. It is found that the waveform is characterized 
by the power-law tail for u > 350M. The power-law index numerically computed is ~ 7.0 and also agrees with that 
analytically derived in [5(j as — (21 + 3) = —7. 

We also computed gravitational waves from a static star. As the stellar models, we adopted models A and C. Since 
the odd-parity matter perturbation is time independent in this case, gravitational waves propagate freely. This implies 



B. Perturbation code 




(6.9) 



(6.8) 
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that with the momentarily static initial perturbation, no gravitational radiation should be emitted. We checked that 
this is indeed the case except for the tiny amount emitted soon after the simulation started, which is due to numerical 
errors associated with the finite differencing and relaxation of the numerical system. 

VII. STELLAR COLLAPSE TO BLACK HOLES AND GRAVITATIONAL WAVES 

We have computed stellar collapse of models B and D to black holes. Since they are marginally stable against 
gravitational collapse, we extracted the internal energy e by 1 % initially to induce the collapse. Here the energy 
extraction is done on the initial null cone. The difference between the energy extractions on the initial null cone and 
on the spacelike hypersurface is very small for nonrelativistic stars but may be significant for highly relativistic stars. 
During the collapse, we adopt the T-law equation of state as 

p=(r a -l)ne. (7.1) 

The simulations were carried out using both the May- White and the Hernandez-Misner schemes. We note that with 
the latter scheme, we can follow the evolution only outside an event horizon. This implies that we cannot find apparent 
horizon and event horizon in the following simulation. We stop the calculation of the stellar collapse and estimate the 
value of the perturbation field on the event horizon using the extrapolation when i? s /2M = 1.01 is satisfied for the 
surface radius R s . We have confirmed that the obtained waveform is not so sensitive to the choice of the criterion. In 
the simulation of the collapse of models B and D, we have observed no evident shock wave until the black hole forms. 

A. Collapse of a supermassive star 

Supermassive stars are quasistable objects of mass 5 x 1O 5 M0 < M < 10 10 Afo and possible direct progenitors 
of supermassive black holes [l4l Irsl] . They quasistationarily contract due to radiative cooling to the onset of radial 
instability 0,E1> resulting in formation of supermassive black holes. The quasistatic evolution of rotating super- 
massive stars was recently investigated by Baumgarte and Shapiro 16]. Taking into account that supermassive stars 
are likely to be rigidly rotating and that the adiabatic index is « 4/3, they clarified that the ratio of the rotational 
energy to the gravitational energy is at most ~ 0.009 at the onset of the gravitational collapse. They also found 
that oblateness of supermassive stars around the central region is very small (see Fig. 1 in justifying that the 

Roche model 15] for rotating supermassive stars is adequate. Shibata and Shapiro [l9| numerically computed the 
collapse of a marginally stable rotating supermassive star to a Kerr black hole in the two-dimensional fully general 
relativistic simulation. They indicated that more than 90% of the stellar mass collapses directly into the black hole in 
the dynamical time scale as in the spherical collapse. These results suggest that if we pay attention only to the inner 
region, the collapse proceeds in a nearly spherical manner. Motivated by this fact, we apply the present perturbation 
analysis to compute gravitational waves emitted during the collapse of supermassive stars. 

1. Spherical collapse 

We adopted model B with T a = (4/3) + 0.00142 to model a supermassive star of mass 1O 6 M . In the numerical 
simulation, we typically take 1000 grid points to cover the supermassive star. For the collapse of model B, with 
the May- White scheme, apparent horizon was located near the center in a late time of collapse, but soon after the 
formation, numerical accuracy deteriorates and as a result computation crashed before the event horizon swallowed all 
the fluid elements. On the other hand, the computation can be continued until the null hypersurface reaches the event 
horizon with the Hernandez-Misner scheme; i.e., the whole region outside the event horizon is computed numerically. 
This clearly indicates the robustness of the Hernandez-Misner scheme. The matching is done on the mass shell within 
which 99.7% of the total mass is enclosed. In the following, we deal with this matching surface as the stellar surface. 

In Fig.|SJa) we display the snapshots of the density profile at selected time steps. In this simulation, the spacetime 
settles down to a static one at u ~ 176760M. The calculation has been stopped at u ~ 176763M or r s ~ 175294M. At 
this moment, we have obtained the redshift (l + z s ) ~ 190.7 at the surface. Since the equation of state of supermassive 
stars is soft, the mass is highly concentrated around the center. In the late stage of the collapse, the increase of the 
central density is accelerated, and hence the collapse proceeds in a runaway manner. (This makes the simulation 
without null formulation technically difficult.) Figure EJb) shows the trajectories of mass shells for u > 176400M. 
Each mass shell asymptotically approaches a constant value greater than 2m because the lapse function decreases 
to zero. This figure shows that the central region collapses earlier, while the outer envelope accretes slowly after 
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the evolution of the central region is almost frozen. Figure [HJc) shows the total mass contained in the high-redshift 
region in which the lapse function is less than 0.1 [or equivalently (1 + z) is greater than 10 where z denotes the 
gravitational redshift]. This region first appears at the center at u ~ 176510M. We find that 80% of the stellar 
mass is swallowed into this high-redshift region within the time interval Am ~ 15M. After the inner region collapses, 
surrounding atmosphere falls into this high-redshift region spending a much longer time ~ 200M. 

2. Gravitational radiation 

For computation of the nonspherical perturbations, a static initial condition for (3 should be given. However, the 
realistic perturbation profile of the initial data set is not clear for the odd-parity perturbation. The purpose of this 
paper is to study the gravitational waveforms during the formation of a black hole qualitatively. Thus, to investigate 
the dependence of the gravitational waveforms on the initial perturbation profile, we gave three kinds of the initial 
data sets as (1) with /3; n ;t =const, (2) with /3; n i t = exp[— (R/R c ) 2 }, and (3) with /3; nit = cxp{— [(R — R s )/R c } 2 }, where 
the scale length of the inhomogeneity, R c , for the matter perturbation is chosen to be R c = i? s /3. For (1), the 
perturbation is uniformly distributed. For (2), the amplitude of the perturbation in the inner region is larger than 
that in the outer region. For (3), the amplitude of the perturbation is the largest near the stellar surface. 

To check the convergence, computations for case (1) were carried out with N = 1000, 500 and 250, where N is the 
number of grid points in the stellar interior. The wavelength of gravitational waves is roughly comparable with the 
stellar size for the early stage and with the stellar size or the size of the formed black hole for the late stage. Since we 
adopt the comoving coordinate, the number of grid points per wavelength is kept almost constant and of the order 
1000 for 1000-zone calculations. 

The waveform, luminosity, and integrated total energy of gravitational waves for I — 2 are plotted in Fig. |5J Since 
the oscillation period of gravitational waves is much shorter than the total time of integration, we show them only 
around the time of major emission. Figure demonstrates that the convergence is achieved well. For example, it is 
found that the totally radiated energy shows the first-order convergence. Although the numerical error in the time 
of the major emission might look very large, it is only due to the extremely long time for the collapse compared to 
the duration of the major emission. It is seen in Fig. Ela) that the waveform is made up of the precursory wave, 
quasinormal ringing, long-time-scale decay, quasinormal ringing again, and power-law tail. Figures EJb) and (c) 
indicate that the primary and secondary contributions to the totally radiated energy come from the first and second 
bursts of quasinormal ringing, respectively. The reason of the existence of this second burst will be explained later. 

We show the waveform, luminosity, and integrated total energy of quadrupole (I = 2) gravitational waves for 
different initial perturbations (l)-(3) in Fig. 1101 In all three cases, the amplitude of the perturbation is normalized 
so that q = 2M. Here, we display the results with N = 1000. 

As seen in Fig. 1101 the gravitational waveform depends strongly on the perturbation profile initially given. For 
case (2), as seen in Figs. HOf c) and (d), gravitational waves of high amplitude are emitted around u ~ 176510M, 
approximately at the same time as the formation of the high-redshift region. The waveform is characterized mainly 
by a quasinormal mode of the formed black hole. The duration of this major emission is roughly ~ 50M, i.e., 
approximately equal to the period and/or the damping time of the quasinormal mode. Indeed, the waveform for 
u > 176520A/ is well fitted by a damped oscillation of the complex frequency 2Mlu w 0.74 + 0.19z, which agrees 
with the theoretically predicted value 2Mui — 0.74734 + 0.17792i [4!| within a few percent error. It is possible 
in principle in the present scheme to observe the modulation of gravitational waveforms due to mass accretion as 
Papadopoulos and Font [Hl| demonstrated for the spherically symmetric Klein-Gordon field on the accreting black 
hole background. However, in the present numerical results, it is difficult to distinguish the effects of the mass increase 
from the gravitational waveforms, possibly because the mass accretion rate is rather high and the duration of the 
accretion is not so long. 

For case (1), as seen in Figs. llOf al and (b), gravitational waves look as the linear combination of a quasinormal 
mode of a black hole and a long-time-scale and nonoscillative component (note that the amplitude of gravitational 
waves does not settle down to zero for a long-time duration ~ 200M after u ~ 176510M). The long-time-scale 
component is produced due both to the perturbation profile initially given and to the nature of the collapse of the 
background spherical star: In the collapse of supermassive stars, the central region collapses first and subsequently, 
the outer envelope gradually falls into the black hole spending a long time duration ~ 200AJ. Since f3- in i t =const, 
the outer envelope retains a considerable fraction of the perturbation for case (1). As a result, the quasinormal 
modes are likely to be continuously excited for the duration ~ 200M during which the matter falls into a black hole. 
However, the quasinormal modes continuously excited should cancel each other due to the phase cancellation effect 
|52|. This suppresses the amplitude of gravitational waves and produces a long-time-scale component, which results 
in the suppression of radiated energy as seen in Fig. HUT g'). 

For case (3), in which the matter perturbation is retained mainly near the stellar surface, as seen in Figs. ITUT c) 
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and (f), the effect of the phase cancellation is more outstanding. In this case, the amplitude of gravitational waves is 
highly suppressed and the amplitude of the quasinormal mode ringing is much smaller than that of the long-time-scale 
component. We note that the time scale of the long-time-scale component is in approximate agreement with the time 
duration in which the accretion of the surrounding envelope continues. 

Figure lTUT h) shows the one-sided power spectral density for the obtained gravitational waveforms for all three cases. 
It is clear that the phase cancellation effects significantly suppress the high-frequency component for cases (1) and (3) 
while the low-frequency component depends not on the initial distribution of matter perturbation but on the initial 
moment of perturbation. 

To explain more clearly why the long-time-scale and nonoscillative component appears, we artificially superimpose 
the waveform <E>(i) obtained for case (2) displayed in Fig. llUf c). which is characterized by quasinormal ringing, with 
the weight factor w(t) as 



dt'w(t')$(t-t'), 



where w(t) is chosen by trial as 



;(t) = 







for < t < t a , 
for t < 0,* a < t, 



(7.2) 



(7.3) 



and we set t c = 2.95M and t a — 236M. The t c and t a correspond to the time scales of the collapse of the inner region 
and of the accretion of the outer envelope, respectively. See Fig. Illf c) for the shape of this function. In the above 
functional form, the first and second terms imitate the effects of inner collapse and subsequent accretion of the outer 
envelope, respectively. In Fig. lllf a1. we display the result for the superimposed waveform. This result is qualitatively 
the same as that obtained for case (1) displayed in Fig. HOf a). We note that the long-time-scale and nonoscillative 
component is produced by the long-term superposition of quasinormal ringing, including a precursory "burst" wave. 
Next we choose another weight factor w(t) defined as 



w(t) 



(t a %) (t.,%) + 



for < t < t a , 
for t < 0, t a < t, 



(7.4) 



and we again set t a — 236M. See Fig. Illf c) for the shape of this function. This functional form imitates the effect of 
accretion of the outer envelope alone. In Fig. lllf b). we display the result for the superimposed waveform. This result 
is qualitatively the same as that obtained for case (3) displayed in Fig. llOf e). It should be again emphasized that the 
superposition not of the pure damped oscillation but of the full waveform for case (2), including the precursory burst 
wave, produces the long-time-scale component. The above consideration confirms our speculation. 

We note that the second quasinormal ringing seen around u ~ 176750M in Figs. HUf a) and (b) for case (1) and in 
Figs. llOf e) and (f) for case (3) are excited by the stellar surface which falls into the black hole finally. Since the phase 
cancellation is not effective at this stage, the quasi-normal ringing dominates the waveform for the latest phase. After 
this ringing, a power-law tail is seen, and the power-law index agrees with the theoretically predicted value —7. 

Finally we note the following point. In this simulation, all the matter is swallowed into a black hole eventually 
because the zeroth-order solution is spherically symmetric. That is the reason why the accretion ends at some finite 
moment. However, in a realistic nonspherical problem, the matter around the stellar surface does not fall into the 
black hole, and would form surrounding disks. This implies that the gravitational waveforms calculated for the latest 
phase u > 176750M will be modified in a realistic simulation. 



B. Collapse of a neutron star 



We performed a simulation for the collapse of a neutron star to a black hole adopting model D. As in the collapse 
of a supermassive star, 1000 grid points are typically taken to cover the neutron star interior. For the collapse of 
model D, both the May- White and Hcrnandez-Misner schemes can evolve the collapse to the formation of the event 
horizon with reasonable accuracy. In the following, the results in the Hernandez-Misner scheme will be described. 
In this simulation, the matching is done on the mass shell within which 96.1% of the total mass is enclosed. In the 
following, we deal with this matching surface as the stellar surface. 

In Fig. E| we show the snapshots of the density profile at selected time steps, trajectory of the mass shells as a 
function of time, and the mass fraction of the fluid elements in the region of strong gravitational field in which the 
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magnitude of the lapse function is less than 0.1. We note that in this simulation, the event horizon is just about to 
form at u ~ 111.5M. The calculation has been stopped at u ~ 111.5M or t s ~ 68.04M. At this moment, we have 
obtained the redshift (1 + z s ) ~ 22.89 at the surface. In contrast with the collapse of the supermassive stars, the 
central density increases only by a factor of 3.0 throughout the collapse. Thus, the evolution of the central density 
does not show the runaway behavior. Also, the collapse proceeds very coherently: Figure IT2T c) indicates that the 
high-redshift region first appears at the center at u ~ 99. 1M and that almost all the fluid elements collapse to the 
high-redshift region in a short time interval Au ~ 6M. 

To study the evolution of the perturbation, we gave the momentarily static initial data sets with matter pertur- 
bations (l)-(3) as in the collapse of a supermassive star. The amplitude of the perturbation is normalized so that 
q = 2M. 

In Fig. 1131 we show the waveform, luminosity and accumulated energy of gravitational waves for case (1). Computa- 
tions were carried out with N = 1000, 500, and 250, and Fig. H3ldemonstrates that convergence is achieved (note that 
results with TV = 1000 and 500 agree well with each other). A small modulation is found just after the beginning of 
computation. This is because the initial data sets which are numerically constructed are not precisely static solutions 
of the finite differencing equation for the evolution of perturbations. The above is confirmed by the fact that this 
modulation becomes smaller as N increases. After this spurious precursor disappears, a black hole quasinormal mode 
is excited and a power-law tail follows |55| . No other outstanding feature is found in the gravitational waveform in 
contrast with the case of the supermassive star collapse. This is likely due to that the collapse proceeds very coherently 
for F a = 2. The first and major gravitational emission of the quasinormal oscillation occurs approximately just after 
the appearance of the high-redshift region. From the numerical result, the quasinormal frequency and the index of 
the power-law tail are calculated as 2Mui = 0.752 + 0.176i and —7.1. These values agree well with the theoretical 
values 2Mw = 0.74734 + 0.17792i and -{21 + 3) = -7 [MIH both within a few percent error. 

Figure ITU shows gravitational waves for different initial matter perturbations (l)-(3). The waveform depends only 
very weakly on the initial perturbation distribution. The reason is that the collapse proceeds in a very coherent 
manner and hence the difference of the perturbation distribution does not generate outstanding differences. There is 
no phase cancellation effects on the neutron star collapse irrespective of the shape of the perturbation as opposed to 
the supermassive star case as is seen in the shape of spectrum. A small time delay of the oscillational phase is seen 
for cases (1) and (3) compared with case (2). This is simply because the matter perturbation enters the black hole 
for case (2) earlier than for cases (1) and (3). 

VIII. SUMMARY 

We have reported a new implementation in linearized Einstein theory. In this code, the Hernandez-Misner scheme 
is adopted to compute a spherically symmetric zeroth-order solution. As a result, we can compute stellar collapse 
to a static black hole until the null hypersurface reaches the event horizon, and the whole region outside the event 
horizon is numerically generated. We emphasize that the collapse of a supermassive star to a black hole proceeds in a 
runaway manner, i.e., the central density grows rapidly although the surrounding atmosphere does not collapse very 
rapidly. It is not technically easy to compute such collapse with a 1+1 numerical scheme. The Hernandez-Misner 
coordinate system is essential to enable thorough computation of the black hole formation. 

We have also proposed a new numerical method to compute gravitational waves in perturbation theory. We divide 
the computational domain into three regions. For computation of the perturbations in the exterior region, we adopt 
double- null coordinates which agree with the characteristic curves of gravitational waves. This choice enables the 
computation of gravitational waves emitted during the entire history of black hole formation. 

To study the qualitative nature of gravitational waves from stellar collapse, we performed simulations for the collapse 
of a supermassive star and a neutron star to black holes. In the gravitational collapse of the neutron star, gravitational 
waveforms are characterized by a black hole quasinormal mode, as demonstrated in the fully general relativistic 
simulations 8J. On the other hand, for gravitational collapse of the supermassive star, the waveform depends strongly 
on the perturbation profile that we give initially. For a centrally concentrated matter perturbation, the waveforms 
are characterized mainly by a black hole quasinormal mode, as in the collapse of neutron stars. However, when the 
matter perturbation is distributed uniformly, the waveform is determined by a linear combination of the black hole 
quasinormal mode and a long-time-scale component which results from the superposition of the quasinormal ringing 
component. Moreover, when the matter perturbation is located around the surface, the long-time-scale component 
dominates the waveform. This is likely due to the less-coherent nature of the collapse of supermassive stars: the 
central part collapses earlier and subsequently the outer envelope accretes on to the central black hole. 

We have determined the choices of initial distribution of perturbation not assuming physically realistic situations 
because of the lack of our knowledge on physically realistic odd-parity perturbation. Therefore, at present, we do 
not claim that the gravitational waves obtained here are realistic in particular for the collapse of supermassive stars. 
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However, our present results strongly suggest that the gravitational radiation is not so sensitive to the initial condition 
for neutron star collapse but is highly sensitive for supermassive star collapse. 

In the formation of intermediate-mass black holes (mass > 2OOM ) formed from quite massive stars, which is 
destabilized by electron-positron pair creation |53j |. collapse would proceed in the same manner as for supermassive 
stars. Thus, gravitational waveforms from the formation of intermediate-mass black holes also depend strongly on 
the state of the precollapse star. 

We have focused on perturbations of odd parity. The dominant modes of gravitational waves are likely to be of even 
parity in most cases, and thus the study of the even-parity perturbations seems to be more important. In a rotating 
stellar collapse, quadrupole deformation that rotating stars retain before collapse will be the source of gravitational 
waves of even parity. The study of such effects on gravitational waves emitted in a black hole formation is now in 
progress. 
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FIG. 1: Spacetime structure of spherically symmetric black hole formation is depicted. Region I denotes the stellar interior, 
region II denotes the intermediate exterior region from which we can emit an ingoing light ray which hits the stellar surface 
before the horizon formation, and region III denotes the far exterior region. The boundary between regions I and II is the 
stellar surface. 
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FIG. 2: Snapshots of the velocity profile in the numerical simulation of homogeneous dust collapse with the initial radius 
R — AM. Crosses and solid curves denote the numerical and exact solutions. The initial data set is put on the initial 
outgoing null cone on which the dust surface is momentarily static. The vertical and horizontal axes are the coordinate velocity 
U = e~^R tU and the circumferential radius R, respectively. The attached labels denote the proper time r s of a comoving 
observer at the stellar surface. All the quantities are shown in units of M = 1. 
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FIG. 3: Gravitational mass as a function of central density of spherical equilibrium stars with polytropic equations of state 
P = Kn T % where (a) T a = (4/3) + 0.00142 and (b) T a = 2. Models A and C are stable while models B and D are marginally 
stable against gravitational collapse. The numerical values are shown in units of c = G = K = 1. 
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FIG. 4: Evolutions of the central density (a) for model A with T a = (4/3) + 0.00142 and (b) for model C with T a = 2. The 
horizontal axis is the time u for an observer at infinity. The internal energy e is initially increased uniformly by 1% both for 
(a) and for (b) from the equilibrium configuration. All the quantities are shown in units of M = 1. 
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FIG. 5: Linear gravitational waveforms for I = 2 in the Minkowski spacetime extracted at r = 5. The initial data set is 
provided on the initial null surface t = r. Three levels of the grid resolution as Ar — Au — Av — 0.001, 0.002, and 0.004 are 
adopted. The solid line denotes the exact solution, (b) is the magnification of the region encompassed by the square in (a). 
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FIG. 6: Total radiated energy of gravitational waves for I = 2 from homogeneous dust collapse for different initial (maximally 
expanded) radii. The initial moment q is normalized so that q = 2M. We set the momentarily static initial data on the 
initial null slice. The results for two kinds of initial distribution for matter perturbation /3i n i t are plotted. For the case of 
,0init oc exp[— (R/Rc) 2 ], the radius R c is chosen to be one-third of the initial radius of the dust ball. The present results are 
compared with those of Cunningham, Price, and Moncrief [2^| (CPM1978), in which the initial data sets were set on the 
spacelike surface of maximum expansion. 
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FIG. 7: (a) Waveform, (b) luminosity, and (c) accumulated energy of gravitational waves for I — 2 radiated from homogeneous 
dust collapse with the initial radius R — 20M. These are estimated at R — 40M and normalized so that q = 2M. The 
horizontal axis is the observer time u. We set the momentarily static initial data on the initial null surface with the matter 
perturbation /?i n i t = const. The solid, long-dashed, and dashed lines denote the results for N = 1000, 500, and 250, respectively, 
where JV is the number of spatial grid points inside the dust ball. The matching is done at the outermost mass shell. All the 
quantities are shown in unit of M = 1. The long-dashed and dashed lines are almost indistinguishable because they lie on top 
of the solid line. 
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FIG. 8: (a) Snapshots of the density profile, (b) the trajectories of mass shells, and (c) the mass fraction contained in a 
high-redshift region are plotted for the gravitational collapse of model B with T a = (4/3) + 0.00142. We extract 1 % of the 
internal energy from the equilibrium configuration to induce the collapse. In (a) the label for each curve denotes the observer 
time u. In (b) the vertical and horizontal axes denote the circumferential radius and observer time, respectively. The labels 
denote the mass fractions enclosed by mass shells. In (c) the high-redshift region is defined as the one in which the lapse 
function a is less than 0.1. We deal with the mass shell which encloses 99.7% of the total mass as the stellar surface. All the 
quantities are shown in units of M = 1. 
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FIG. 9: (a) Waveform, (b) luminosity and (c) accumulated energy of gravitational waves for I — 2 radiated from the collapse of 
model B with T a = (4/3) +0.00142. Gravitational waves are extracted at R = 2000M. The initial matter perturbation is chosen 
so that ,0i n it = const. The momentarily static initial perturbation is given. The amplitude of the perturbation is normalized so 
that q — 2M. The solid, long-dashed, and dashed lines denote the results for N = 1000, 500, and 250, respectively, where N 
is the number of spatial grid points inside the star. The matching is done on the mass shell which encloses 99.7% of the total 
mass. All the quantities are shown in units of M = 1. 
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FIG. 10: (a), (c), and (e) Waveform, (b), (d), and (f) luminosity, (g) accumulated energy, and (h) the one-sided power spectral 
density of gravitational waves for I = 2 radiated from the collapse of model B with T a = (4/3) + 0.00142. Gravitational waves are 
extracted at R = 2000M. The momentarily static initial perturbation is given. The amplitude of the perturbation is normalized 
so that q = 2M. In (g) and (h), the solid, long-dashed, and dashed lines denote the results for /3i n it = const, exp[— (R/R c ) 2 ], 
and exp{ — [(R — R B )/R C ] 2 }, respectively, where /3i n it is the initial distribution of the renormalized matter perturbation and _R C 
is chosen to be R c = R s /3. For waveforms and luminosities, we display the results for each case in each figure for clarity. Note 
that the vertical axis of (g) is in logarithmic scale. The matching is done on the mass shell which encloses 99.7% of the total 
mass. All the quantities are shown in units of M — 1. 




FIG. 11: Waveforms generated by superimposing the waveform obtained for case (2). The shape of weight factors are displayed 
in (c). (a) and (b) show the waveforms for the weight factor labelled "Collapse+ Accretion" and for that labelled "Accretion." 
The waveforms are normalized so that the maximum amplitude is unity. See text for details. 
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FIG. 12: Same as Fig. [8] but for the collapse of model D with T a = 2. We deal with the mass shell which encloses 96.1% of 
the total mass as the stellar surface. 
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FIG. 13: Same as Fig. [5] but for the collapse of model D with T a = 2. Gravitational waves are extracted at R = 100A/. The 
matching is done on the mass shell which encloses 96.1% of the total mass. All the quantities are shown in units of M = 1. 
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FIG. 14: (a) Waveform, (b) luminosity, (c) accumulated energy, and (d) one-sided power spectral density of gravitational 
waves for I = 2 radiated from the collapse of model D with T a = 2. Gravitational waves are extracted at R = 100M. The 
momentarily static initial perturbation is given. The amplitude of the perturbation is normalized so that q = 2M. The solid, 
long-dashed and dashed lines denote the results for /3i n i t = const, exp[— (R/R c ) 2 ], and exp{— [(R — R B )/R C ] 2 }, respectively, 
where /3i n i t is the initial distribution of the renormalized matter perturbation and 7i c is chosen to be R c = R s /3. The matching 
is done on the mass shell which encloses 96.1% of the total mass. All the quantities are shown in units of M — 1. 



